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We show how to construct relevant famihes of matrix product operators in one and higher di- 
mensions. Those form the building blocks for the numerical simulation methods based on matrix 
product states and projected entangled pair states. In particular, we construct translational in- 
variant matrix product operators suitable for time evolution, and show how such descriptions are 
QQ ■ possible for Hamiltonians with long-range interactions. We show how those tools can be exploited 

' for constructing new algorithms for simulating quantum spin systems. 
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The study of strongly correlated quantum systems is currently receiving a lot of attention. To a large extent, this 
is due to the formidable progress that has been made in creating such systems under controlled laboratory conditions 
, such as in optical lattices and ion traps. From the theoretical point of view, major new insights have been obtained into 
■ characterizing the nature of the wavefunctions associated to those strongly correlated systems. The concept of matrix 
' product states and its generalizations plays a central role in those new insights, as it provides a sound foundation 
f~| and justification for the success of numerical renormalization group methods and especially of DMRG Those 
Oh insights have led to the development of new algorithms for simulating quantum spin systems; most notable are the 

. algorithms for simulating time evolution [1, 0, S 13 and the ones generalizing DMRG to higher dimensions . 
^ ' In this work, we are concerned with the efficient construction of so-called matrix product operators (MPO), the basic 
^ building blocks for those novel algorithms. MPO were introduced in the paper 0, Q and form the operator analogue 
of matrix product states. We will show how to construct translational invariant MPO's in 1 and 2 dimensions that 
I— I, approximate real or imaginary time evolution; in contrast to the TEBD/DMRG algorithms the translational 

symmetry is not broken in the Trotter step. This generalizes the constructions reported in [l^. Second, we construct 
MPO descriptions for general Hamiltonians with decaying long-range interactions. This is very interesting in the light 
of simulating quantum spin systems with long-range interactions. 
I Similar work for constructing MPO representations of Hamiltonians has independently been reported in 
0^ ■ Reference ^ gives a very nice presentation of matrix product operators from the point of view of DMRG, and also 
contains results on how to write spin chain Hamiltonians using MPO. Reference explores the connection between 
matrix product operators and Markov processes in depth, and also contains some results on generalizations to higher 
[ dimensions. In reference an algorithm is devised to simulate quantum spin chains with long-range interactions in 
00 I the thermodynamic limit; it also contains similar results as reported here on the approximation of power law decaying 
interactions by sums of exponentials. 
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MATRIX PRODUCT OPERATOR DESCRIPTIONS OF EXPONENTIALS 

A. Construction 



Let us first start with a simple example: suppose we want to simulate real/imaginary time evolution under the 
Ising Hamiltonian in transverse field 



<ij> i 

where only nearest neighbour interactions are considered; both the one- and two-dimensional case will be considered. 
As usual, this evolution will be approximated using a Trotter expansion, but we want to do this is such a way that 
the translational invariance in not broken. Therefore, we split the Hamiltonian in two parts Ti. = Hz + where Hz 
contains all terms with operators and H^ the ones with . Obviously, all terms in Hz commute, and therefore 
Oz — ex-p{eHz) can be calculated exactly. As we will show, Oz has a very simple and elegant MPO description, and 
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of course Ox has a trivial MPO description as it is a product of strictly local operators. Time evolution can now be 
described within the formalism of matrix product states by evolving the MPS under the action of the MPO 0^0^. 

Let us next show how the MPO of Oz can be constructed. First, observe that 



exp(eZ(g)Z) = cosh(e)/ (8) / + sinh(e)Z (g) Z 



Z(g) z 



Here we used the notation Z'^ = I , = = Z and defined the vectors Bi. Let us now consider the translational 
invariant 1-D case of N spins with periodic boundary conditions 



exp {^^Z.iZ,+i 



= JJexp(eZiZj+i) 

i 

E {{BlB,,){BlB,J...iBlB,,)) Z^Z^ ® Z\^ Z 



E Tr {Bj,BlB,,BlBj,...B,,B]^) Z\-+'- ® Z 
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All A;2... 



E Tr {C^^C^K..C^'') ® Z. 



In the third step, we made use of the cyclic property of the trace, and in the fourth step, we made a change of variables 
fci = i\ © ji where binary arithmetic is assumed. We have therefore proven that exp (e ZiZi^x) has a very efficient 
matrix product description with the matrices given by 



' ' I smh(e) / 

„i _ tVT _ ( v^sinh(e) cosh(e) \ 

C - - 1, Vsinh(6)cosh(e) j 



A big advantage of this precise MPO formulation is that is is symmetric; the spectral properties of the associated 
transfer operator are hence well behaved, which is important if used in algorithms with periodic boundary conditions 

In two dimensions, we can repeat exactly the same argument and obtain the PEPS description of the operator 



exp e E = E F{C''\C='\ ...)Z^' ® Zf ® ... 

y <ij> j xxx^... 

with tensors 

Clp^,8= E BMB,{P)Bk{l)Bi{5). 

i-\-j^k-\-l—x 
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Here Bi{a) means the a component of the vector Bi, x € {0,1} and the sum is taken over i,j,k,l = : 1 with 
the condition that i+j + k + l= x in binary arithmetic. This proves that the PEPS description of the operator 

exp (^£j2<ij> ^i^j^ bond dimension 2. Note that no approximations were made and as such this statement is 
vahd for any value of e. In particular, this gives the MPO description for the classical Ising partition function; its free 
energy can therefore be calculate by contracting the tensor network consisting of tensors C*'. 

The previous analysis can trivially be generalized to the case of any Hamiltonian that is a sum of commuting terms: 
for this class of Hamiltonians, exp(eff) has a very simple matrix product operator description. As this holds for any 
e, it also holds for all thermal states, and by taking e — > — oo it is proven that all ground states of such Hamiltonians 
have exact MPO descriptions that can easily be constructed. Notable examples of this is the toric code state of Kitaev 
and the family of string net states [3, [l^ . 

From numerical considerations, it is useful if the matrices/tensors occurring in the MPO description are real and 
symmetric. There are some tricks of how to achieve this. Consider for example the Heisenberg antiferromagnetic 
Hamiltonian 



Ti-Heis — ^ {XiXj + YiYj + ZiZj) . 

The operator exp {—f3Ti.Heis) can be decomposed in Trotter steps consisting of H^, Hy, H^, and every Trotter term 
involves operators of the form exp(— eiJ^;). As we saw in the previous section, the associated matrices involve terms 
like sinh(e), which becomes complex when e > 0. What we can do however is a change of basis on every second site 
(this obviously only works for bipartite lattices), where we rotate the spins with the unitary operator Y = ay] this 
maps — ^ —X2n, Y2n ^2n j ■^2n — > — ■^2n- On thc Icvcl of the Hamiltonian, this flips the sign of the and 
interactions, for which the associated operators expl+eHx) have indeed real and symmetric MPO descriptions. The 
problem seems to remain however with the operator exp(— ei?y). This can however easily be cured by defining the 
real antisymmetric matrix Y — iY for which Hy = —Hy when we replace all operators Y by Y. Next, exjp{+eHy) can 
again be expressed as a MPO; however, we have to be careful as Y.Y — ~I as opposed to +/. Looking back at the 
derivation of the MPO for the Ising case, we can easily see that this sign can be absorbed into C, and we can express 



\ i / feife2... 



exp ( -e y^YY+i ) = > ^ Tr ( C'^C^.. ) Y^'' ® Y^ 

as a MPO with matrices 

cosh(e) 



^sinh(e) cosh(e) 

■\/sinh(e) cosh(e) 



Of course the same can be done in two dimensions. Here we get 



where the sum in the power y— 1 is not in binary arithmetic. This clearly leads to a real translational 

invariant MPO parametrization. 



B. Algorithms 

It is now obvious how to turn those MPO-descriptions to our advantage for constructing new algorithms for the 
simulation of quantum spin chains. 

Let us first consider the case of imaginary time evolution, where the goal is to evolve a state in imaginary time such 
as to simulate a thermal (finite /3) or ground state (Ji oo). Obviously, we will use the Trotterization described in 
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the previous section. The big advantage there is that the translational invariance is never broken, and furthermore 
that the matrices involved in the MPS description of the MPO are real and symmetric. In particular, that means 
that, if we start with a translational invariant MPS with real symmetric MPS description, then it will stay like that 
during the whole course of the evolution. This has a dramatic effect on the numerical conditioning and stability of 
the algorithm. 

The algorithm for time evolution is now as follows: given a translational invariant MPS with matrices {A^} with 
bond dimension D and MPO with matrices {B^}, {X*} of dimension we want to find a way of representing cutting 
the bond dimension of the MPS {C*} given by 

C ="^A^ (g) B''{i\X^\j) 
jk 

in an optimal way. This can easily be done as follows: calculate the leading eigenvector x of the transfer operator 
E = C** <8i C* (note that E is symmetric and as such this is a very well conditioned problem). Rewriting a; as a 
DD' X DD' positive semidefinite matrix, we can easily calculate its singular value decomposition x — UT,U^ . We now 
define the projector/isomctry P as the rectangular matrix consisting of the first D rows of U, and act with this P on 
the matrices C^ The updates matrices A* are therefore obtained by = P^C^P which is obviously still symmetric 
and real. Clearly, all those steps have to be done in such a way as to exploit the sparse nature of the problem, such 
as done in DMRG, which leads to a complexity that scales like D^. Also, if the eigenvalues that are thrown away are 
not small enough, we can always increase the bond dimension. 

The big advantage of this procedure is that it is extremely well conditioned and very efficient to implement. This 
allows to go to very large bond dimensions. Notably, as compared to the TEBD algorithms, we do not have to 
take inverses at any time (because the gauge degrees of freedom are trivial as they consist of unitary matrices), and 
furthermore it works equally well if the MPO is very far from the identity operator (this is important in the context 
of PEPS algorithms). 

The same ideas can of course be used in the case of real time evolution. In that case, the matrices involved become 
complex symmetric, and it might be benefitial to apply some gauge conditions to optimize the stability. This can be 
done as follows: given x, we want to find the complex (symplectic) matrix Q, Q.Q^ = 1 such that the condition number 
(i.e. smallest divided by largest singular value) of the matrix QxQ^ is as large as possible. This optimization problem 
can be solved recursively as follows: calculate the singular value decomposition oi x — vsv'', choose the generator 
Gk = —G^ as Gk = lui{vDv\y — viv\), and make the substitution x — s- exp(— ieG)a; exp(ieG) for small enough e, and 
repeat this until convergence. Convergence is equivalent to the derivative of the condition number being zero. The 
final gauge transform to be implemented is the product of all infinitesimal transformations Q = J| exp(ieG'fe). Note 
again that all of this becomes trivial in the case of real symmetric matrices (such as occurring in imaginary time 
evolution): in that case the Q cannot change the condition number as they are unitary. 

We have tested those new algorithms on the critical Ising and Heisenberg spin chain models, and obtained results 
that are consistent with what we expected. In particular, for the Heisenberg antiferromagnetic spin chain, we obtain 
a precision of (i?D=64 — Eexact) / Eexact — 2.83 * 10^^ with very modest calculations. In the case of the critical Ising 
in transverse field, we get (i?_D=64 - E act) / E exact = 1, 10 * 10^9. 

The algorithms for the 2-D analogue will be discussed elsewhere [l7| . 



II. MATRIX PRODUCT OPERATOR DESCRIPTIONS OF HAMILTONIANS WITH LONG-RANGE 

INTERACTIONS 

A. Construction 

Let us next investigate how to represent Hamiltonians with long-range interactions of the form 

with f{i — j) some decaying function. The first question to ask is whether it is still possible to find an exact MPO 
description of O = exp(eH). It can easily be seen that this is not possible if the function f{x) does not vanish at 
some finite distance: otherwise, the action of O on a MPS could increase the Schmidt number over any cut with an 
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arbitrary large amount, and hence no finite MPO description is possible. This is the reason why the transfer matrix 
approach in classical 1-D spin systems breaks down for such long-range interactions. 

So let's be less ambitious and try to find a MPO description of the Hamiltonian itself. This is interesting for several 
reasons: first, this is useful in constructing algorithms for time evolution using iterative methods like Lanczos, and 
second, it allows to calculate quantities like {tp\H'^\ilj) efficiently. 

As a start, let us consider a general 1-D spin 1 /2 Hamiltonian with nearest neighbour interactions. If the Hamiltonian 
is translational and reflection invariant, then there always exists a basis such that the Hamiltonian can be written as 

a.i j 

where O can be any one-qubit operator. Similarly to the construction of MPS descriptions of the W-state [l^l, a 
MPO can be constructed to represent this Ti. by making use of nilpotent matrices: 



iii2--- 

Xq — I Xi = <Tx X2 — <Jy X3 ~ az Xi — O 

VI = |0) Vr = |4) 

Bo = |0)(0| + |4)(4| 

Bi = |0)(l|+/ii|l)(4| 52 = |0)(2|+Ai2|2)(4| B3 = |0)(3|-fAi3|3)(4| 

B4 = |0)(4| 

The simplest way of deriving this is to think about the Hamiltonian a Markov process with 5 possible symbols 
(remember that MPS can be constructed using Markov processes), such that a symbol Xi, X2, X-}, is always followed 
by itself and then all zeros Xq, and X^ by all zeros. As such, one can easily prove that I? = 5 is optimal in this case 
because this is the operator Schmidt number of the Hamiltonian when splitting it into two pieces. Note that if only 
Ising interactions would have been considered, then _D = 3 would have been sufhcient and we could have chosen 

i3o = |0)(0| + |2)(2| Bi = |0)(l|+Mi|l)(2|. 

Note that there is no need for B2,Bs, B4 in that case. 

It is obvious how to generalize this description to the case of higher dimensional systems and to the case of 
exponentially decaying interactions. Let us first look at the case of exponentially decaying interactions. By adding 
diagonal terms to Bq 

Bo = |0)(0| + A,|l)(l| + A,|2)(2| + A.|3)(3| + |4)(4| 
we can immediately check that that the corresponding Hamiltonian / MPO is given by 

a,i<j j 

which is a spin chain with exponentially decaying interactions. 

Unfortunately, it is impossible to get exact MPO descriptions when the interactions are decaying following a power 
law. However, inverse polynomials can pretty well be approximated by sums of exponentials (this is the reason why 
DMRG is able to reproduce the correlations in critical models pretty well). Hamiltonians with power law decay 
of correlations can therefore be well approximated by sums of MPO's, which is itself a MPO. Actually, very few 
exponentials are needed to get a good approximation, even at large distances. The problem of finding the optimal 
weights and exponents for such an approximation problem for a general function /(fc), i.e. 



N 

min^|/(fc)-^x,A,'=|, 



' k=l 



is not completely trivial. In the appendix, we present a simple method that solves this optimization problem for 
general f{k) and a given number of exponentials n and a number of sites N > n (the method works for any functions, 
and returns complex exponents in the case of oscillating functions as should be). If we choose power law decay with 
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cube power 3, N=1000 and n=10 then the above cost function is 10 ^ (the maximal difference between the function 
and the approximated one is 5.10~*). This maximal difference falls to 3.10"^ for power 2 and 3.10^^ for power 1. 

In conclusion, we found the exact MPO description for Hamiltonians with exponentially decaying interactions. 
Hamiltonians with power law decay can be approximated very well using sums of such MPO. The matrix product 
operators obtained for the description of Hamiltonians are of a very different form than the ones obtained by taking 
the exponential. The main difference is that the corresponding transfer matrices will always contain a Jordan block 
structure, and one has to be careful in dealing with such situations when considering the thermodynamic limit. 

Let us now turn to the 2-dimensional case. Let us again first consider the square lattice with only nearest neighbour 
interactions. There is a very simple way of writing down a PEPS description that achieves the task: first, consider 
the MPS 

\w) = j2m^\o)..\omk\o)...\o)^. 

k 

which is the equal superposition of having one spin up and all other ones down over all sites. Note that this MPS has 
bond dimension 2, and can therefore trivially be represented as a PEPS with bond dimension 2. The idea is that this 
excitation specifies where to put an interaction. Let us next consider the tensors 

Bl<.,p,,,s = |0) (00001 

4aA7,5 = |1)(00|((01| + (10|) + |0)((01| + (10|)(00| 
where we assume that the indices a, /3 are the left respectively top indices, and the associated operators 



= J 

= Z 

It can readily be seen that we get the Ising Hamiltonian if we act with the state on the fifth index of the 
tensors: the state puts one index i equal to one, and the other terms are such that an interaction to the right 
and below it will be created. The total bond dimension of the corresponding PEPS (including the is therefore 
4. 

Decaying interactions between one spin and all other ones can however be obtained in a much more elegant way; as 
we will show, it is even possible to model power law decay of interactions exactly. The idea is as follows: the critical 
classical Ising model in 2 dimensions has power law decay of correlations. Consider the quantum state 



\i,p) = exp -/? Z,Z, (1+))^ 



where |+) stands for the superposition |0) + This is obviously a PEPS, as it is obtained by acting with a MPO 
(see earlier) on a product state. The partition function of the Ising model at temperature /3 is obtained by calculating 
the overlap 

((+l)^^IV'(/3)), 

and correlation functions between two spins are obtained by replacing the corresponding (+| at the left side of this 
expression by (— | = (0| — (1|. Instead of the \yV) state in the previous example, we will use a state that is the 
equal superposition of two excited spins as opposed to one. The MPS description of has bond dimension 3 and 
is similar to the ones derived for the 1-D Hamiltonians with exponential decay where we put the parameter A = 1 
(i.e.Bo = J^)- This gives us all the necessary ingredients to construct the MPO description of the Hamiltonian: 

n = ((a;'^|(a;'^|...(a;'''|) (|V/3)|W^"))X'' 0X'^...(g)X'« 

iii2...ijv 

= |0)|+), ki) = |l)|-) X^=I X' = Z 

Here the vectors act on two qubits, one on the corresponding qubit of \'ipp) and the other one on The 
\W") state enforces that exactly two operators X"^ will be nontrivial, and jf/'/s) gives the right weight to the associated 
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interaction. As those are products of PEPS, the result is a PEPS with bond dimension 2x3 = 6. This is am amazing 
result: as opposed to the 1-D case, there is an exact PEPS description for Hamiltonians with 2-body interactions that 
decay as the with 1^ = 1 the critical exponent of the Ising model. 

Obviously, this construction can be repeated for any classical spin model, and hence many different exponents can 
in principle be taken. It is as yet an open question how to engineer the PEPS such as to get a specific exponent, 
although very good approximations can again be obtained by making use of sums of exponentials. 



It is obvious how to make use of all this in algorithms for simulating quantum spin chains. First of all, it is clear 
how to extend the variational MPS method described in [l^ to the present case. For this, we have to consider finite 
systems with open boundary conditions and matrix product states that have site-dependent matrices in their MPS 
description. The optimization {iplH]^) can then be done using the alternating least squares method described in p^ . 
As expected, numerical tests showed very good convergence properties. 

The problem of time evolution is a little bit more challenging, as we cannot use the Trotterization tricks. However, 
Krylov-based methods can of course be used (see [l^ for a review), and are the method of choice here. 

These finite dimensional algorithms can not readily be generalized to the infinite case however. A sensible way for 
determining ground state energies in that limit would be to first use the finite dimensional algorithm just described, 
and then use a brute-force gradient-based optimization method for optimizing the infinite case. For this we need to 
be able to calculate expectation values {tp\0\ip) in the thermodynamic limit, i.e. when 1-0) is an infinite MPS with 
matrices {A^} and O the MPO description of a Hamiltonian with exponentially decaying interactions. The idea is 
to consider a family of MPO On whose support is limited to N sites (i.e. the Hamiltonian does only act on N 
sites), calculate the energy with respect to the infinite MPS (this energy will scale linearly in N), and the take the 
thermodynamic limit to calculate the energy per site. It holds that 



where \L) = \xi)\vi),\R) = with \xi), \xr) the left/right eigenvectors of the transfer matrix Eq = J2i^i ® 

the vectors w/,fr are the ones used in the MPO description of the Hamiltonian, and 



The eigenstructure of Eh is nontrivial because it has Jordan blocks. In the present case, the only relevant blocks 
are of size 2 (larger blocks would lead to an energy that scales superlinearly with the size of the support of H, which 
cannot be). Using the notation used before, one sees that the left/right eigenvector corresponding to the largest 
eigenvalue in magnitude do is given by {qi \ = (a;;|(0| / \qr) = \xr)\0). The generalized eigenvectors can now be found 
by solving the equation {Eh — d^I) \qr) = |gr-), [Eh — d^I) = {qi\. We next define the matrices Qi = {\qi), I?/)), 
Qr = {\<lr)^ \<lr)) and the 2x2 matrix Q = (QfQr) ■ In the limit of large N, it holds that 



B. Algorithms 



Eh =^^A,(E) Bj (g) Ak{k\Xj\i). 



ijk 




As QT 




Q'^\qi), the expectation value ('0|Ojv|'0) is given by 




in the limit of large TV. The energy per site is therefore given by 



Ql2 



^/{qi\qr) 



l/{qi\{EH-doir^\qr). 



In a similar vein, it is possible to calculate expectation values of the operator {ip\{H — A)^]?/"). This is relevant 
because it gives an exact bound on how far a given MPS is from an exact eigenvector of H. As we have a squared 
term, Jordan blocks of dimension 3 will be encountered. As before, we define 
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Eh^ ="^Ai(g>Bj(giBk(E> Ai{l\xlXj\i). 

ijkl 

The relevant right eigenvector is again of the form \qr) = |a;r)|0)|0) and we can find the eigenstructure of the 
associated Jordan block as follows: start by calculating \qr) as we did in the previous section using the operator Eh 
instead of Ejj2; next define \(]') = sym{\qr)\0)) where the symbol 'sym' means symmetrization with respect to the 
part of the state acting on the Hamiltonian part of the MPO (the antisymmetrized wavefunction turns out to be an 
irrelevant eigenvector of Eh2 with eigenvalue rfo); finally solve the linear set of equations {Eh^ — dg/) |g"r) = l^r)- 
The relevant eigenstructure is now given by the matrix Q,. = {\qr) , Wr) ' W r)) ■ Similarly, one can calculate Qi and 
Q = {QjQr)~^- The final expectation value is then given by 



( 1, N, N{N - l)/2 )QyOj 

The matrix Q therefore contains all the relevant information about the energies and their scaling when N — > oo. 

The energy can now easily be optimized with a brute force gradient-base optimization routine. 

Concerning the 2-dimensional MPO representing Hamiltonians, it turns out that they are very valuable for speeding 
up actual calculations done by the PEPS method: the calculation of the expectation value of the Hamiltonian with 
respect to to a given PEPS can be calculated in one run using this idea, and we don't have to calculate the expectation 
value for every terms individually anymore. 



III. CONCLUSION 



In conclusion, we constructed several examples of interesting families of matrix product operators in 1 and 2 
dimensions. Those descriptions turn out to be very valuable for constructing stable and scalable algorithms for 
simulating quantum spin systems, in one and two dimensions. 



IV. APPENDIX 



In this appendix, we show how to solve the problem of approximating any function f{k) as a sum of exponentials 
for k = 1..N: 



N 



i=l 



min^|/(A:)-^x,A,^|. 

First, construct the rectangular N — n + lx n matrix 

/ /(I) /(2) /(3) 
/(2) /(3) ... 
/(3) 



f{n) \ 



\/(iV-n + l) 
/ A? 

A? 



A^ 
A^ 



f{N-l) 

m-i) f{N) I 

^° ^ /xi 

X2 



iV-ra \N-n 



\xr" A^ 



\ 



A^-y 



V 



A? 
A? 


A} . 
A^ . 


• An 

• A^ 


A° 




/ 



w 
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Note that is a Vandermonde matrix. We observe that F and W span the same space (note that A'' is typically 
much larger than n), such that there exists a, nxn matrix Q s.t. FQ ~ W. Define Fi as the rectangular matrix which 
consists of the first N — n rows of F and F2 as the one with the last N — n rows. Due to the Vandermonde structure of 
W , it must be approximately true that FiQK ~ F2Q with A the diagonal matrix containing the exponents. Therefore, 
A ~ Q~^F^^ F2Q {F^^ denotes the pseudoinverse of Ui): the exponents {A;} hence correspond to the eigenvalues of 
the matrix F^^ F2 which can be calculated very easily. 

This method can be made more robust by making use of the so-called QR-decomposition. This can be done by 
first calculating the (economical) QR decomposition of = UV and by defining Ui as the rectangular matrix which 
consists of the first N — n rows and n columns of U and U2 as the one with the last N — n rows: there must again 
exist a Q such that UQ ~W . The exponents A can therefore easily be calculated as the eigenvalues of the matrix 
U^^U2- The advantage of using the QR-decomposition is that the pseudoinverse of Ui is much better conditioned 
than of Fi . 

Once those exponents are found, a simple least squares algorithm can be used to find the corresponding weights 
{xi}. It happens that this method is very efficient and reliable, even when oscillating functions are involved. A similar 
method is known in the field of signal processing under the name of Hankel singular value decomposition. 
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